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A thermohydrodynamic lattice-BGK model for the ideal gas was derived by Alexander et 
al. in 1993, and generalized by McNamara et al. in the same year. In these works, particular 
forms for the equilibrium distribution function and the transport coefficients were posited 
and shown to work, thereby establishing the sufficiency of the model. In this paper, we 
rederive the model from a minimal set of assumptions, and thereby show that the forms 
assumed for the shear and bulk viscosities are also necessary, but that the form assumed for 
the thermal conductivity is not. We derive the most general form allowable for the thermal 
■ conductivity, and the concomitant generalization of the equilibrium distribution. In this way, 

we show that it is possible to achieve variable (albeit density-dependent) Prandtl number 
even within a single-relaxation-time lattice-BGK model. We accomplish this by demanding 
j-jjfy analyticity of the third moments and traces of the fourth moments of the equilibrium distri- 

bution function. The method of derivation demonstrates that certain undesirable features of 
the model - such as the unphysical dependence of the viscosity coefficients on temperature - 
cannot be corrected within the scope of lattice-BGK models with constant relaxation time. 

o 
o 

> 

A thermohydrodynamic lattice BGK model for the ideal gas was derived by Alexander et al. in 1993 (!]]. 
In the same year, McNamara and coworkers Q generalized this model to include different relaxation times 
for the stress and the energy, thereby allowing adjustment of the Prandtl number. These analyses assumed 
particular forms for the equilibrium distribution function, took the transport coefficients to be linear in the 
internal energy / and independent of the density p, and carried out a Chapman-Enskog analysis to show 
that these choices yielded the desired thermohydrodynamic equations. In this way, the sufficiency of the 
thermohydrodynamic lattice-BGK model was established. In this paper, we rederive this model from the 
demand that the moments of the equilibrium distribution be analytic. We show that the forms assumed for 
the shear and bulk viscosities in earlier works are necessary as well as sufficient, but that the form assumed 
for the thermal conductivity k is not. We derive the most general possible form for k, as well as the alteration 
that this incurs in the form of the equilibrium distribution function, and we show that this enables some 
adjustment of Prandtl number even within a single-relaxation-time lattice-BGK model. 

We demonstrate this by carrying out an "inverse Chapman-Enskog analysis" ; that is, we assume only that 
the first three moments and the trace of the fourth moment of the equilibrium distribution function exist 
and are analytic, and we demand that the resulting macroscopic equations be identical to the usual ones of 
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thcrmohydrodynamics. We then work backwards to find the requirements imposed upon the equilibrium and 
the transport coefficients. As shall be shown, this process determines the derivatives of the third moments 
and the traces of the fourth moments of the equilibrium distribution function with respect to the conserved 
densities. Analyticity of the moments then requires that their mixed second derivatives be symmetric with 
respect to the order of differentiation. 

It is sometimes glibly stated that this "inverse Chapman-Enskog" approach to deriving lattice Boltzmann 
models is capable of dealing with essentially any hydrodynamic equations whatsoever. This work shows 
that this is definitely not the case. In fact, as we shall see, even the problem for the ideal gas is hugely 
overdetermined. Such overdetermination is a serious problem with lattice BGK models in spite of the 
freedom that one has in choosing different lattices and different sets of particle velocities. Indeed, within 
the context of a lattice BGK model for the usual equations of viscous thermohydrodynamics, with constant 
relaxation parameter, this paper will demonstrate that no choice of lattice or of velocity set will yield 
viscosities that vary as the square root of temperature, as is usually desired. This follows from nothing more 
than the requirement that the moments of the equilibrium distribution be analytic. It may be possible to 
circumvent this restriction by other generalizations of the BGK methodology - for example, by introducing 
a particular dependence of the relaxation parameter on the hydrodynamic densities - but this has yet to be 
demonstrated. 



II. THE CHAPMAN-ENSKOG ANALYSIS 



As usual, we assume that the mass, momentum and kinetic energy densities are given by the first three 
moments of the discrete-velocity distribution function, 

i i 

P = E C *^ =E c */i 0) 

i i 

F - V °< f - \^ C i f (0) 

i i 

where f\ is the equilibrium distribution. The kinetic equation can be written 

/ i (x + ec i ,i + eA*) = / < (x,t)+n i , (1) 

where we have denoted the lattice vectors by Cj , and the time step by At which is equal to unity in natural 
lattice units. We have also introduced the gradient expansion parameter e which will aid in the Chapman- 
Enskog analysis, and the collision operator £li which is required to conserve mass, momentum and kinetic 
energy 

E°i = ° 

i 

E = 

i 

E f = °- 

i 

The assertion that the energy is purely kinetic will lead us inevitably to an ideal gas equation of state. 
Throughout this work, we assume that the collision operator is of BGK form, 

0* = \ {fi 0) ~ fi) • ( 2 ) 

where r is a characteristic collisional relaxation time. As noted above, this model was generalized in subse- 
quent work by McNamara et al. to include more than one relaxation time and thereby allow for variable 
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Prandtl number. In this work, we adhere to the original form with a single relaxation time. We shall show 
that some variation of Prandtl number is possible even in this case. We emphasize, however, that there is 
nothing preventing the application of our methods to the more general case of multiple relaxation times, and 
that such application may uncover still more general lattice-BGK models for that case as well. 

To carry out the Chapman-Enskog analysis, we first solve Eq. (0) for the full distribution function 
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where Di = dt + v, • V and v, = Ci/At, and where Eq. (Q) was used to get from the first line to the second. 
Next, we establish multiple time scales by ordering Di as follows: 



Di X'" ' D 
fc=i 



fc-i n( fc ) 
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where 

Here ti is the fastest time scale, £2 is slower, etc. This results in 
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and we can expand this in e to obtain 



fc=0 



where to C(e 2 ), 
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or, in expanded form 



f (i) 



f(2) _ 



-r^+v.-V)/! 05 
-t [(St 



v.-V)-(r--](a tl +v t .Vr 



f (0) 



(3) 
(4) 



III. THE FORM OF THE EQUILIBRIUM DISTRIBUTION 



To proceed, we adopt notation for the first five completely symmetric tensor moments of the equilibrium 
distribution function. Since the local equilibrium is entirely determined by the conserved quantities, all of 
these moments must be functions of the three conserved densities, thus 
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Q a p(p,P,£) = ^c ia c ifj fl 0) 
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where Greek letters denote spatial indices. We also adopt notation for the traces of these tensors, 

Q aa {p, p,£) = 2£ 
R a <xp{p,P,£) = 2Kp(p,p,£) 
S aa /3j(p,P,£) = 2Lp 7 (p,p,£). 

We have not shown and shall not show in this paper how to find a functional form for an equilibrium 
distribution with these moments. Suffice it to say that the need to be able to do this places well understood 
constraints on the lattice symmetry group |3). This is why, for example, a triangular lattice is often used in 
two dimensions, and a face-centered cubic lattice is often used in three. Henceforth, we shall simply assume 
that the lattice is such that moments of the distribution transform like isotropic tensors to rank four. 

IV. THE FIRST-ORDER SOLUTION 

Taking moments of Eq. (^), we get the conservation equations 

d tlP + V-p = 
d tl P + V-Q = 

d tl £ + V K = 0. (5) 

These equations tell us how the conserved densities vary on the fastest (ti) time scale. Since all of the 
moments are assumed to be functionally dependent on the conserved densities, we also know how they vary 
on the fastest time scale. 



V. THE SECOND-ORDER SOLUTION 

Taking moments of Eq. (Q), we get the conservation equations 

d t2 p + V ■ p = (r - I) V • [dip + 2ft, V • p + VV : Q] 



2, 
V 

V 

Using the first-order solutions, Eqs. (^), to simplify the above, and converting to index notation, we have 



dt 2 P + V ■ Q = ( r - 1 ) V • [3 t 2 lP + 20 tl V Q + VV : R] 
d t2 £ + V • K = (t - i J V • [d\£ + 2d H V K + VV : L. 



4 



d t2 p + Vf3 P = O (6) 

dt 2 Pa + VpQpa = - ^ V/3 (d tl Qff a + V 7 i? 7/3ct ) (7) 

9 t2 £ + VpKp =( T ~l) v /3 (fti^/J + V 7 L 7/3 ) • (8) 

To proceed, we make use of the assumed functional dependence of Q and K on the conserved densities, 
p, p and 8. Using the chain rule and Eqs. (|j|), we have 

dtxQpa = Q/3a,pdt 1 P + Q (ja,p- l dt 1 p 1 + QfSa^dtiS 

= — <2/3a,pV 7 p 7 — Q/3a,p 7 V 'sQs-f ~ Q/3a,sV 7 -K'.y 

= — Q,3a,pV 7 p 7 — Qf3a,p y (Q8j,p^SP + Qs-y,p e ^5P(, + Qs-f,£^s8) 

-Q/3a,e (-^7,pV 7 p + K^^V^ + K^^V^E) 
= — (Qfla^-yQs^.p + Q/3a,£Ks,p) VsP 

— (Q/3a,pSs^ + Q/3a,p^Qsy,p ( + Q/3a,£^<5,p 5 ) V ' gpg 

— (Q/3a,p- 1 QS'y,£ + Q/3a,£-K5,f) V $8 , 

and likewise 

d tl Kp = - (K ^Qs ltP + Kp t£ K SiP ) V s p 

— (Kf3, p 8s£ + K / 3 tPi Qs 7 ,p i + Kp t£ Ks yPi ) S7sP^ 

— (Kp,p^Qs-y,s + KpfiKsfi) V58 , 

where the commas denote differentiation [J]. Inserting these into Eqs. (^) and (ft), we get 
d t2 p a + VpQ Pa = ( T ~l) v /3 {vZs^sp + Vf aSt V m + Vf a5 v s 8 
8 t2 8 + V P K = (t - ~\ Vp (V%V SP + V^VsPt + Vf s V 5 8 ) , (9) 
where we have defined the transport coefficients 



^PaS = R$l3a,p ~ Q /3a,p-,Q 6-y ,p ~ Q/3a,£Ks, 
■^3at5£ = R$f3a,p£ Q f3a,p3s£ Qpa,p-,Qs~f,p£ Q /3a, £ ^<5,pj 
^$aS = R&Pa,£ ~ Qj3a,p^Q&^,£ ~ QfiafiKf,^ 

T> £ p P s = Lsp,p — Kp lPy Qsj, P — Kj3 i sKs tP 
^p% = Lsp, Pi — Kp tP 6g£ ~ K l 3^ Pi Qs 7 ,p i — Kp t sKs, P( 

Dps = Ls[3,£ - Kp jP ^Qs lt £ - Kp.gKgx. (10) 

These are the most general equations that we can obtain using the lattice Boltzmann method that we have 
described. Note that Q and K are determined by the desired momentum and energy fluxes. Eqs. ([To]) then 
give us differential equations for the moments R and L in terms of the conserved densities. 



x Thus, for example, Kp tP = dKp/dp. 
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VI. DESIRED THERMOHYDRODYNAMIC EQUATIONS 



The equations of thermohydrodynamics can be written in conservative form as follows 

= dtp + V • (pa) 

= d t (pu) + V ■ (puu - P) 

= d t {pe) + V • (pue u P - kVt) 

where p is the density, u is the hydrodynamic velocity, e is the total energy per unit mass, and 

u 2 

is the internal energy per unit mass. Also, we have defined 



(11) 
(12) 
(13) 



P = —PI + A (V • u) 1 + /i Vu + (Vu) J 

where P is the pressure. 

Comparing the conserved densities in the above equations with those used in the last subsection, we see 
that the meaning of p is unchanged and that we must identify 

p = pu 

£ = pe. 

We also define the internal energy density, 



I = pi = £ 



V 



and assume henceforth that the pressure P(p, I) is a function only of the mass density and the internal 
energy density. Recasting Eqs. (11) through (|13|) in terms of p, p and £ is then straightforward, and after 
some algebra we find 



dtPa + V/3 

d t £ + V/3 



where we have defined 
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V/3 
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A- 



Pa_ 

P 



£ Vip Spsp^VsPt + kSpsVgS 



(14) 



and 



A 
k 



X 



(r-l/2)p 
P 

(r-l/2)p 
k 

(r-l/2)p 



That is, we have incorporated the factor of r — 1/2 into the definitions of the barred transport coefficients 
to facilitate comparison of Eqs. ( |l4] ) and (|9|). Henceforth, we assume that A, /2 and k are, like the pressure 
P, functions of the mass density p and internal energy density / only. 



G 



We can now identify the required momentum and energy fluxes, 



n P0Pa , pA 

W8a = r rOpa 

P 

K (i = ^{£ + P). (15) 
P 

Thus, the second moments, Q a f3, and the trace of the third moments, R aa f3 — 2Kp, are completely deter- 
mined in terms of the conserved densities. To proceed, we have to show that third and (traces of the) fourth 
moments with the required properties exist. Necessary and sufficient conditions for this are discussed in the 
next subsection. 



VII. CONSISTENCY CONDITIONS 

To determine the tensors R and L, we compare Eqs. ( O ) and (||) to make the identification 

(16) 

(17) 
(18) 

/>> VaVi T kSfiS ( f P 2 \ /i n\ 

v 0s = -—r^m ~~^~\ £ ~j) ( 9) 

&m = —Km - —fa (20) 



T>PP 

U paS 


P 


T)PP 




U /3aS 


= 



m = hs 0S . (2i) 



Hence, from Eq. (nfj) we have 



Rs(3a,p = Q/3a,p T Q(57,p + Qf3a,£Ks,p —A-a/38^ 

Rsf3a, Pi = Q{3a,pSs£ + Q 0a,p~,Q 5^ ,p ( + Q/3a,£Ks iP( + h. a j35£ 

Rspa,£ = Q/3a,p^Qsj,£ + Qpa^KsX, (22) 



and 



L S0,p = ^-!3,p^QSj,p + i^l3X K S.,p ~2~ AasPSt, — \ t — 

Ls/3,pt: = Kp. p 5s£ + K^ jP ^Qs JlPs + Kfj^Ks, P( + —A a psz ~~^ a<s 

Lspx = Kp. p ^Q Sl:£ 4- Kp^Ksx + k5p S . (23) 

Since Eqs. (|l5|) give the components of Q and K in terms of conserved densities, it follows that the entire 
right-hand sides of Eqs. ( p2] ) and ( p3| ) are known in terms of the conserved densities. Eqs. (|2^) and ( |23| ) 
then give the partial derivatives of R and L, respectively, with respect to the three conserved densities. 
Thus, analyticity of R and L requires that all mixed second derivatives be symmetric; in addition, we must 
demand that the trace of R be equal to 2K, as required. The determination of the second derivatives is 
straightforward in principle but very tedious, being somewhat simplified by the assumption that P, A, p, and 
k depend only on p and /. The results are presented in Appendix [A]. A glance at this appendix shows that, 
even accounting for the isotropy and symmetry under index interchange of Rsp a and Lgp and of the second 
derivatives with respect to the momentum density, there are dozens of potentially independent requirements 
to be satisfied. It is therefore something of a minor miracle that all of these requirements can be reduced to 
the following six scalar conditions 
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A + P P Pj 
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I + P 



P,pJ + P, 



p,p 



I + P 



P \ P 
I + P 



PjJ + P,pJ 



P.i + P.i 



= A - fi + k + Ik j + pk yf . 



(24) 
(25) 
(26) 
(27) 

(28) 
(29) 



in any number of dimensions. The first pair of these, Eqs. (24) and (25), must be integrated for /2; analyticity 
then imposes a consistency condition on the mixed second derivatives of ft, but it is seen to be satisfied 
identically, and these equations then integrate to yield 



/' 



P 
P ' 



(30) 



Likewise, the second pair, Eqs. (|26|) and ( |27| ) for A, are automatically consistent with analytic A, and integrate 
to yield 



(31) 



The fifth equation, Eq. fl28|), is then satisfied identically. This leaves Eq. (|29|) which is a partial differential 
equation for the thermal conductivity; we shall return to this equation in a moment. 
With the aid of Eqs. ( ^0|) and (|3l|), Eqs. (Q) can be integrated to get 

o rx , x , x \ P , PBPpPa 

ttgpa = (Of3aPS + OSaP/3 + OspPu) 1 5 ■ 

P P 

For D — 3, this result may be compared to Equation (8d) in the reference by McNamara We must now 
demand that the trace of this quantity be equal to twice the heat flux. We find 



Rgpp — 



P v 2 
P P 



Ps, 



where D = 8 aa is the spatial dimension. Equating this with twice the heat flux, as given by Eq. (|l^), we 
arrive at 

which is the ideal gas law 0. From this, Eq. (^) immediately yields 

21 



P = 



Dp' 



and Eq. (ml) yields 



2 To see this in a more familiar form, note that the internal energy per unit mass is / = ^pkuT /m, so that P — nksT 
where n = p/m is the number density. 
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AI 



Returning to the thermal conductivity, it is straightforward to establish, e.g. using the method of charac- 
teristics, that Eq. (E9I) has the general solution 



2(D 



D 2 p 



I + f 



where /' is (the derivative of) an arbitrary function of its argument. Thus, we find that the model does 
have some flexibility in the choice of thermal conductivity. Using this result, it is now straightforward to 
integrate Eqs. (p3[) for the trace of the fourth moment, and the result is 



L 



50 



2{D + 2) 
Dp 



I 



Dp 2 



2{D + 2)p 



I + Pf 



Sg/3 



D + A 
D + 2 



Dp 2 



2(D + A)p 



PSP/3 
P 



For D = 3 and / = this may be compared to Equation (8c) in the paper by McNamara [MJ. To the best 
of our knowledge, however, this form for L has never been written down for general D, or with the / term 
that allows some control of the thermal conductivity. 

The "unbarred" forms of the transport coefficients are then, 



p 

A = 
k = 



2{D 



D 2 



I+f 



We note that p and A are independent of p, while k depends on it only via the arbitrary function /. We see 
that the relation 

, 2 

holds; when D = 3 this is approximately true for many real gases. The linear temperature dependence of A 
and p are, however, unrealistic; for real gases these go as \fl. Finally, the Prandtl number is 



Pr^ = 



D 



D 



1 



f'(I/p) 



For / = 0, this reduces to the constant D /(D + 2) which is indeed equal to 1/2 for the D — 2 example treated 
by Alexander et al. Q|, and 3/5 for the D = 3 example treated by McNamara |0J (when the relaxation times 
in the latter paper are set equal). More generally, however, the function / may be used to exercise some 
control over the Prandtl number. 



VIII. CONCLUSIONS 



We have presented a first-principles derivation of the lattice-BGK model of the ideal gas. In the process, 
we have established the necessity of the functional forms used for the shear and bulk viscosities, and have 
derived a somewhat more general form than is usually assumed for the thermal conductivity. We have also 
shown that the thermohydrodynamic model of the ideal gas is grossly overspecified, so that its existence is 
largely fortuitous. We have also shown that the unphysical dependence of the viscosities on temperature is 
an essential feature of the lattice-BGK model with constant relaxation time. 

The cornerstone of our demonstration was nothing more than the demand of analyticity of the third 
moment and the trace of the fourth moment of the equilibrium distribution function. The main point of this 
paper is that it is essential to pay close attention to this criterion when using the "inverse Chapman-Enskog" 
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method of derivation of lattice Boltzmann models. It is a very fundamental and restrictive requirement 
which cannot be circumvented by different choices of lattice or of velocity set, or other superficial details. 
As noted earlier, allowing the relaxation parameter to depend on the hydrodynamic densities may be a way 
to recover some flexibility in this regard, but this possibility remains completely unexplored. 

One of the principal outstanding problems in lattice Boltzmann research is the development of a thermo- 
hydrodynamically consistent model that includes a soft interaction potential . The ultimate solution to 
this problem will surely involve a radical departure from current lattice Boltzmann models - for example, by 
inclusion of some information about the two-particle distribution function. However it is solved, we expect 
that the form of the equilibrium distribution function involved will not be as easy to intuit as it is for the 
ideal gas. In that event, we expect that an approach similar to the one that was adpoted in this paper will 
be useful. 
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APPENDIX A: MIXED SECOND DERIVATIVES OF THE MOMENTS 



The mixed second derivatives of the third moments and the trace of the fourth moments of the equilibrium 
distribution function can be calculated as described in Secton 7. The results are 



Rsfic 



,p,Pa 



Rsfja, 



($6f3S aii + SsaSfe) P,p + 



P - Pf 



X + P p Pi , fI + P\ 



+ 



+- 



5&fj [ PaPt - y<W ) + Ssa [ ppPt - 



Pj 



3a \ PSPi - yC 



P 



PpJ + P,P,P 
- 

2 ' 
I + P 



Pj - 



P.JJ + P,pJ 



Rspa,p,£ — Rs0a,£,p 



PS 



5pc 



p V p J 



P.JJ + Pj,p 



~- (SaSPp + 5(38Pa) \ Pj 



Pj_ 

P 



Rs/3a, Pi , Pn - R8f3a,p, 1 , Pi = +~ [(SspS ar) + Sp v 5s a ) P^ ~ {Ssa$P£ + Ssp5 a ^) p v ] ^/J,j -^j 

I'.l.l + PpJ 



P\ 



I + P 



Rst3a, Pi ,£ — Rsp a .,£, Pi — +^8^f)c 



A,, + 3 + (^)p iJiJ + PiJi , 



Lsp.p.p^ — Lsf3,p i ,p — + 



1 



P- 
P 2 , 

-FT PJ 



pj 



+ 



2p 
h^pp ( 
P 2 I 



k + Ik j + pk tP — fj, — P/U l(C 

_Pj 
P 

p - (^r~) {Pj + pPp l) + pX ' p ~ pp 



+PiP P + pP, P ,p] 



p 



p 2 I 



-p PiPjPk 



A + 
A 



- p 2 i 

X.I + — + 



Pj + pp„ 



I + P 



2p 



PjJ + PpJ 



P, 



7T~ I Pj 



9 



Pj fl + p\ 



+ P P j + PJ-^y 



Sp,p,£ - Lsp,£,p = 



_&sp_ ( 

P I 

PiPj 



k + Ik j + pip - (^y^j Pj -P P +j (pj ~y)} 
f , P J , ( I + P \ 



Pjj + P, P j + Pj 



S5^p v - dg v 6p£ 



A- p + 



I + P\ 



9 



Pj + P, t 



11 



+ 



PS {Sp-qPZ ~ SfcPr,) 



+ 



P/3 {$6 V P£ ~ 8 S £P V ) 



P\ 



1 
P 



Lsf3,p ( ,£ - Lsf3,£, Ps — +- {SspPz + Sp^ps) [ H,i 



+ 



_ p] 

A j H — — + 

p 



A / H — — + 
P 

P.j 

P 

fI + P\ 



I + P 



r.i.i + p p j 



/'/./ + p p j 



12 



